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Abstract 

The origin of ferroelectricity in KH2PO4 (KDP) is studied by first-principles 
electronic structure calculations. In the low-temperature phase, the collective 
off-center ordering of the protons is accompanied by an electronic charge der- 
ealization from the near and localization at the far oxygen within the O-H- • O 
bonds. Electrostatic forces, then, push the K + ions towards off-center positions, 
and induce a macroscopic polarization. The analysis of the correlation between 
different geometrical and electronic quantities, in connection with experimental 
data, supports the idea that the role of tunnelling in isotopic effects is irrelevant. 
Instead, geometrical quantum effects appear to play a central role. 



1 Introduction 

Among the hydrogen-bonded ferroelectrics, the compound KH2PO4 (KDP) has been 
the most intensively studied over the second half of the last century. 0, ||] However, 
a detailed knowledge about its ferroelectric phase transition and some related phe- 
nomena is still lacking. In particular, its transition temperature (T c « 122K), shows 
an increase of about 107K upon deuteration. The origin of this huge isotopic effect 
is still controversial. In the sixties, Blinc and Zecks proposed a quantum tunnelling 
model to explain it, which focuses purely on mass-dependent effects. ||, f§] Improve- 
ments of this model include coupling to vibrations of the K-PO4 complex. 0] Experi- 
mental studies over the past fifteen years, especially high-resolution neutron difraction 
measurements, |5|, |[ 0, || provided increasing evidence that the geometrical modifi- 



cation of the hydrogen bonds upon deuteration, the well-known Ubbelohde effect, [[TC 
is intimately connected with the mechanism of the phase transition. McMahon et al. 
have shown that the distance 5 between the two collective equilibrium positions of 
the protons is remarkably correlated with T C .|J They proposed that 5 might be the 
only relevant parameter for T c and, in addition, no concluding evidence of tunnelling 
was established. It was then suggested the necessity of a revised theory for this class 
of transitions. Recently, a phenomenological model was proposed, which couples, via 
geometrical details of the H-bond, the proton coordinate with a nonlinear dynamics 
of the K-PO4 system, where the spontaneous polarization actually develops. [^] This 



model, which includes both tunnelling and geometrical isotope effects, also allows to 
explain the T c behaviour upon deuteration. 

Summarizing, it appears well-established that, in this class of systems, proton or- 
dering and geometrical and electronic properties of the host cage are connected in a 
non trivial way. ||12| So far, to the best of our knowledge, no ab initio calculations have 
been reported for KDP. This type of approach has the advantage of allowing a confi- 
dent parameter-free analysis of the microscopic electronic and structural changes upon 
proton ordering that drive the ferroelectric transition. To this end, we have performed 
Density Functional (DFT) electronic structure calculations [13|, 14 for KDP, within 
a gradient corrected approach to exchange and correlation, both in the tetragonal 
paraelectric phase (forcing the protons to be in the center of the H-bonds) and in the 
orthorhombic ferroelectric phase, where the protons order cooperatively at off-center 
positions in the H-bond. We thus analyzed in detail the structural parameters and 
the microscopic electronic density redistributions associated to the proton ordering 
phenomenon, which eventually provoke the ferroelectric instability. 

This work is organized as follows: in Section 2 we describe the methods and the 
details of the calculations. Section 3 is devoted to the analysis of the geometry results. 
In Section 4, we study and discuss the polarization mechanism produced by the off- 
centering of the hydrogen atoms along the H-bonds, and its relation to the ferroelectric 
instability. Finally, in Section 5 we elaborate our conclusions. 



2 Ab initio methods 

We have performed two different types of ab initio calculations. First, we have used a 
recently developed, very efficient localized basis (LB) set approach (SIESTA). This is 
a fully self-consistent, pseudopotential DFT method that employs a basis set of pseu- 
doatomic orbitals (PAO).[[TR |T(| We have chosen a double- zeta basis set with polariza- 



tion functions (DZP). The exchange-correlation energy terms were computed using the 
Perdew-Burke-Ernzerhof gradient-corrected functional. ||17|| Angular-dependent, norm- 



conserving Troullier-Martins pseudopotentials |1| were employed to represent the in- 



teraction between ionic cores and valence electrons. We have also included nonlinear 
core corrections (NLCC) for a proper description of the K ion, due to the large overlap 
of core and valence charge densities. We used a 125 Ry energy cutoff for the grid com- 



putation of numerical integrals |U| |15| , and a value of E c = 50 meV for the orbital 
confinement energy, which gives total energies and geometries of sufficient accuracy 
for this system. Due to the approximate character of these confined basis orbitals, 
we have checked the overall accuracy by performing standard pseudopotential plane 
wave (PW) calculations. The parallel version of the Car-Parrinello SISSA FPMD code 



19| , and the same set of pseudopotentials, exchange-correlation functional, and NLCC 
were used. The energy cutoff for the plane wave expansion was set to 150 Ry. In both 
approaches, the electronic Brillouin zone sampling reduced to the T-point, which is 
a reasonable approximation for such a wide band gap insulator. We have, however, 
checked that this sampling is sufficient for the large supercells considered. 

In the paraelectric phase, the KDP structure (space group I42d) has a body- 
centered tetragonal primitive lattice. For our calculations, we used instead the conven- 
tional face-centered tetragonal cell (space group F4d2),f2Cj which includes 64 atoms (8 



Table 1: Comparison of the ab-initio (LB and PW) lattice parameters and internal coor- 
dinates with experimental data. The notation is the same used in the experimental works 
referred. The different cases considered in the calculations are explained in section 2. We 
defined Azk as the displacement of the K atom along z, respect to its centered position 
between P planes. Distances are in A and angles in degrees. 
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formula units). With this choice, the transformation to the orthorhombic structure of 
the ferroelectric phase is simply described by a change in the ratio of the basal plane 
lattice parameters, in addition to changes in the internal structure parameters. 

3 Geometry results 

In order to obtain the equilibrium configurations relevant to ferroelectricity in KDP, 
we performed a series of computational experiments, which are described below. 

First we considered a static view of the paraelectric phase by fixing a tetragonal cell 
with parameters at their experimental values at 293 K |, and relaxed the atomic 
positions under the constraint that the H atoms remain in the middle of the H-bonds. 
In this way we found the transition state (TS) which connects the two oppositely 
polarized realizations of the ferroelectric phase. In a second time we still kept fixed 
the tetragonal lattice, and allowed the hydrogen atoms to relax from their centered 
positions. When the unstable TS state was slightly perturbed, two H atoms approach 
each PO4 tetrahedron. The two PO4 oxygens which were approached by the hydrogens 
(in the following named 02) lie on one side of the (xy)-plane containing the P ion. 
The other two oxygens (01) lie on the other side. The geometry relaxation is carried 
out also for all the other atoms in the unit cell. As a consequence of proton ordering, 
also the K ions move away from the centered position between P planes, thus leading 
to a ferroelectric phase with tetragonal lattice, which we denote FP t . Finally, we 



Table 2: Mulliken populations for each atom in the FP t and TS cases for the LB calculation. 
The last line shows the charge difference Aq = g(FP^) — q(TS) between both cases. Units 
are electrons (e) for the first two rows, and e/1000 for the last one. 





01 02 P K H 


FP t 
TS 

Ag/1000 
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6.201 6.201 4.755 0.208 1.116 
+82 -58 -8 -3 -17 



fully optimized atomic positions and cell parameters. This leads to an orthorhombic 
ground state struture GS, which corresponds with the actual experimental ferroelectric 
structure. 

The resulting structural parameters are shown in Table 1, together with experi- 
mental data for the paraelectric tetragonal and the ferroelectric orthorhombic phases. 
An excellent overall agreement is observed for the TS and FP t configurations obtained 
with the PW method, as compared to experimental data. The good quality of the LB 
method can be assessed by comparison with the PW results in columns 4 and 5. In 
the FP( calculation, the off-center ordering of the hydrogen atoms along the H-bonds 
leads to inequivalent 01 and 02 oxygens. Consequently, (a) there is a distortion of 
the PO4 tetrahedra (see the different P-Ol and P-02 distances on Table 1) and (b) 
the K ions abandon the centered positions between the phosphate planes (see Az K 
in Table 1), thus inducing a spontaneous polarization which leads to ferroelectricity. 
These effects will be explained in the next Section. We also verified that, when the 
hydrogen atoms are constrained to remain in the middle of the H-bonds, then the 
centered positions are stable for the K ions. This provides a strong evidence that the 
origin of ferroelectricity is in the off-centering of the protons in the H-bonds. 

It is worth mentioning that, in this fixed tetragonal FP t calculation, a substantial 
anisotropic stress appears in the xy-plane. Note that the hydrogen off-centering in the 
H-bond network, which lies nearly on the xy-pla.ne, leads to an asymmetry of the x and 
y directions. Further relaxation, obtained by optimizing also the cell parameters in 
order to minimize the stress with respect to ambient conditions, releases the tetragonal 
symmetry in favor of the orthorhombic structure shown in column 3 of Table 1. 

4 Polarization mechanism for the ferroelectric in- 
stability 

In order to understand the origin of the ferroelectric instability, we have focused our 
analysis in the electronic charge redistributions that happen in going from the centered 
(TS) to the off-centered hydrogen case (FP t ). To this purpose we compute the orbital 
Mulliken populations, which are shown on Table 2 for the different atoms in both 
cases. It is known that Mulliken populations depend strongly on the choice of the 
basis set. Differences, however, are much less sensitive, and thus meaningful. In going 
from TS to FP t , it can be clearly observed an increase of the charge localized around 
01, the main contribution (« 70%) being provided by the 02 charge decrease. The 



Table 3: Overlap Populations between atoms in units of e/1000. Only the most significant 
Aq = g(FPt) — <?(TS) values are shown. 
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analysis of bond and self atom-orbitals overlap populations, [^TJ which are shown in 
Table 3, also helps to discriminate the origin of this charge redistribution. 

The trends observed in Tables 2 and 3 are confirmed by the charge density difference 
Ap(r) = p FPt (r) — pTs( r ), which is plotted in the planes determined by the atoms P- 
02-H (Fig. la) and P-Ol- • H (Fig. lb). Analysing both, Table 3 and Fig. la, we 
observe a significant enhancement of the 01-01 overlap population, accompanied by a 
smaller increment of the same sign in the Ol-P orbitals. This happens at the expenses 
of the population of the overlap orbitals in the 01- • -H bond. This behaviour around 
the 01 atom indicates that, as the hydrogen atoms move away from 01, the strongly 
covalent Ol-H bond weakens, its charge being redistributed. A substantial amount 
of this charge localizes around 01, while part of it localizes in the Ol-P orbitals, 
consistently with the shortening of the Ol-P bonds reported in Table 1. The contrary 
occurs in the vicinity of the 02 atom: as H moves towards 02 along the 01- • -H-02 line, 
the 02-H distance decreases, thus increasing the degree of covalency and producing 
a charge localization along the corresponding 02-H bond. This is indicated by the 
bond overlap population in Table 3, and also by the contours in Fig. lb. Finally, the 
02-02 and 02-P orbitals loose a significant amount of charge. We name this process 
of charge redistribution, produced by the collective hydrogen off-centering along the 
H-bond, the polarization mechanism (PM). The PM can be described in a simplified 
way as a charge accumulation and localization around 01 at the expenses of a charge 
depletion and derealization in the vicinity of 02. 

This charge accumulation around the 01 atom can also be viewed as a negative 
charge defect, which electrostatically attracts the positive K + ion towards the PO4 
cage, along the z-axis. This electrostatic force is counterbalanced by repulsion forces 
from orbital overlaps, which stabilize a new equilibrium position of the K + ion, as 
can also be seen in Table 1. Our calculation then supports the idea that proton 
off-centering and ferroelectricity are very correlated phenomena. 

Cell relaxation - from FP t to GS - also affects the internal structure (Table 1). 
As H approaches 02, 5 increases, the 01- • -H bond weakens and the H-02 bond be- 
comes stronger. As a consequence, a larger amount of charge localizes around 01, 
and delocalizes away from 02. This behavior is related to the PM and will be useful 
to establish a connection between S and the charge difference Ago = q(01) — q{02). 
Moreover, we observe a striking correlation between Aqo, the distortion of the tetra- 
hedra, A(P-O) = d(P — 02) — d(P — 01), and the displacement of the K atoms, Azk- 
These trends are shown in the theoretical curves in Fig. 2 for two different values of 
S, obtained from the cell-unrelaxed FP 4 and the GS calculations. 

The crucial experimental observation that the effect of isotopic substitution can 
be almost exactly reproduced by applying pressure [^|, [|, supports the idea that the 
main effect of deuteration is to modify the internal geometry. In this sense, our 
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Figure 1: Differential charge density contours Ap(r) in the planes containing the following 
atoms: (a) P-Ol- • H, (b) P-02-H. Labels 01 and 02 denote the positions of the respective 
nuclei, positioned at (0,0). Labels 02-P and Ol-P indicate the position of the center of the 
corresponding bonds. The same convention is used for the 02-H and 01- • H bonds. Positive 
(negative) contours are in solid (dashed) lines. The thickest lines represent an absolute value 
of 2.96 x 10 -3 eA -3 . The thinner lines are obtained by successively halving this value, down 
to 3.70 x 10~ 4 eA" 3 . 



calculations for FP t and GS, corresponding to two different pressures, can be used 
to analyse isotopic effects. It is known from experiment that H(D) ordering, Azk, 
A(P-O), and spontaneous polarization P s , show very similar trends just below T c , 
in the ferroelectric phase. |22 Moreover, these trends are quite well explained by the 
present PM (left curves in Fig. 2). Therefore, it is natural to attempt to establish a 
connection between Aq Q and structural data measured in KDP and in its deuterated 
analog, DKDP. This is done in Fig. 2, by including the experimental values of Azk, 
A(P-O) and P s for KDP and deuterated DKDP. The comparison between theoretical 
and experimental curves is very favorable, especially for what concerns the slopes of 
the curves. [f23| 

In connection with the proposed PM, the above indicates that the enhancement 
of the saturated polarization due to deuteration, experimentally observed below T c , 



22j| can be explained merely by the increase in the charge unbalance between 01 



and 02 (Ago), which in turn is induced by an increase in 5. This picture is in line 



with the ideas developed in Ref . ||24|| , where the experimental behaviour of P s and the 
inverse dielectric constant were rationalized in terms of a proton-dipole interaction 
model. In fact, at temperatures well below T c , the protons are strongly bound to 
the neighboring PO4 dipoles, in a stable minimum of an asymmetric potential. In 
this asymmetric configuration, the protons localize into the deepest well because the 
energy needed to be promoted to the excited state is too high, thus ruling out the 
possibility of tunnelling. [53] In this case, deuteration can only induce the indirect 
effect of enhancing 5 [p5[| , i.e. the so-called geometric effect. Preliminary ab initio 
calculations indicate, even more strikingly, that there is actually no double minimum 
structure in the potential energy surface for single protons. Appart from supporting 



this scenario, the current ab initio-dehved PM establishes the importance of the short- 
range proton - P04-dipole interaction for the transition, as was also speculated in ||. 
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Figure 2: The variaton of A(P-O) (circles), Azk (triangles), Aqo (diamonds) and saturated 
P s (squares) with the H(D) displacement 5. Left curves: LB theoretical results for cases FP^ 
and GS (open symbols) and right curves: experimental values for KDP with 5 e ^ p ~ 0.34A, 
and for DKDP with b e £ v ~ 0.45A, from Ref . p2| , (full symbols), both in the ferroelectric 
phases at T C -20K and T C -10K, respectively. Experimental 5 values considered are estimations 
of these magnitudes at T c , as they cannot be measured in the ferroelectric phase. [ 22 1 The 
lines are guides to the eye only. Arbitrary units are used in the vertical axis. 



5 Conclusions 

We presented the results of ab initio DFT calculations for the H-bonded ferroelec- 
tric material KDP, which reproduce very well the structure of the tetragonal and 
orthorhombic (ferroelectric) phases in KDP. On this basis we identified the micro- 
scopic polarization mechanism that leads to the ferroelectric instability. This arises 
from an electronic charge redistribution between the oxygen atoms involved in the 
H-bonds, occurring when the protons move away from the center of the bonds. We 
indirectly showed that the interplay between proton ordering and ferroelectric polar- 
ization qualitatively explains the enhancement of the saturated P s upon deuteration. 
In addition we provided, for the first time, fully ab initio evidence supporting the idea 
that isotope effects reduce, at least at low temperatures in the ferroelectric phase, to 
the indirect geometric effect of enhancing 5, as speculated in recent works ||. 
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